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Abstract 


This research explores the viability of using a navigation system that relies on 
measurements of the magnetic anomaly field as an alternative to GPS navigation. 
Previous research has been conducted on developing a navigation system using the 
intensity of the Earth's magnetic anomaly field as an alternative signal. This research 
focuses on using vector and tensor measurements, as opposed to scalar measurements 
of the anomaly field, as a means of obtaining accurate position and orientation solu- 
tions. 

This paper presents two navigation systems. The first uses an Extended Kalman 
Filter (EKF) with vector measurements of the magnetic anomaly field to aid an 
inertial navigation system (INS), while the second uses tensor measurements. 

Simulations examine the performance of both navigation systems in sixteen sce- 
narios. The parameters evaluated in the simulations include the position and velocity 
of the trajectory, whether vector or tensor measurements are used, the quality of the 
INS paired with the filter, and the map resolution. Simulations demonstrate that 
the tensor measurement filter paired with a navigation-grade INS performed best out 
of the sixteen test cases. For a one-hour ship trajectory, the navigation system was 
able to demonstrate 35.94 m DRMS error when paired with a navigation-grade INS. 
The same navigation system was able to obtain navigation accuracies of 38.10 m 
DRMS when paired with a 10X-grade INS for a 25 hour ship trajectory with a lower 


resolution magnetic field map due to the depth of the ocean. 
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NAVIGATION USING VECTOR AND TENSOR MEASUREMENTS OF THE 
EARTH’S MAGNETIC ANOMALY FIELD 


I. Introduction 


This thesis focuses on using a simulation trade study to determine the viability 
of using vector or tensor measurements of the Earth's magnetic anomaly field to 
navigate. Two navigation systems are presented and their ability to obtain an accurate 


estimate of vehicle position and attitude throughout a trajectory is evaluated. 


1.1 Problem Statement 


The Global Positioning System (GPS) has become a staple in civilian and military 
navigation systems. The United States has become highly reliant upon the GPS signal 
given its high accuracy. However, GPS is not without disadvantages. GPS is not 
always available in environments that have features that obstruct the signal, such as 
in heavily wooded areas, or in cities where the user is surrounded by large structures. 
The GPS signal is also susceptible to jamming or spoofing. In these cases, especially 
during military operations, the need for high-accuracy navigation still exists, and 
alternative signals need to be used to achieve this. 

Similar to the GPS signal, the magnetic anomaly field of the Earth is a globally 
available signal. The Earth’s magnetic field is always available and is not as sus- 
ceptible to jamming or spoofing as GPS. Maps of the anomaly field signal are also 
available in varying resolutions over the entire globe [31]. Additionally, the Earth’s 
magnetic anomaly field is a vector field that provides more than one signal to use 


for navigation. This research explores the navigation accuracy that is possible when 


1 


using the directional components of the magnetic anomaly vector field as well as its 


spatial gradients (tensor) to navigate. 


1.2 Research Objectives 


Two navigation systems are presented in this research. One system aids an In- 
ertial Navigation System (INS) by matching measurements of the three directional 
components of the magnetic anomaly field to existing maps of the directional compo- 
nents of the field. The second aids the INS by matching measurements of the spatial 
gradients, or tensor, of the magnetic anomaly field to existing anomaly field tensor 
component maps. 

Both navigation systems are tested in this research through simulation, using 
actual ship and airplane trajectories. Trade-space analysis was done to determine 
how changing different variables affects the accuracy of the navigation systems. The 
variables explored include the position and velocity of the truth trajectory, the types 
of measurements used in map-matching, the quality of INS that was used, and the 
map resolution. The results of this analysis will shed light on the viability of the 


proposed navigation systems as alternatives to GPS navigation. 


1.3 Overview 


This thesis is organized as follows. Chapter II provides a basic introduction to 
concepts used throughout the document, such as modeling and mapping the Earth’s 
geomagnetic field, geomagnetic measurements and the instruments used to collect 
them and a brief background on the prevalent navigation filters used for magnetic 
navigation. Chapter II also provides an overview of previous and related work done 
in the field of magnetic navigation. 


Chapter III outlines the two specific navigation systems used in this research. 


It begins by defining the full dynamics and measurement models of the navigation 
system that uses vector measurements of the Earth's magnetic field. It then defines 
the full dynamics and measuremnet models of the navigation system that uses ten- 
sor measurements. Chapter III also presents details on the simulation framework 
developed, such as the specific trajectories used, the magnetic maps used, and the 
measurements and INS data that was simulated. 

Chapter IV presents the results of the simulation trade study. The navigation 
accuracies of the two different navigation systems were compared for each trajectory 
and notable results were highlighted. 

Chapter V concludes the thesis by giving an overview of the simulation trade 
study results. Chapter V outlines how navigation system performance was affected 
by changing the simulation parameters. Chapter V also includes possible areas for 
future work to improve navigation accuracies of the filters presented in this document 


and to determine their viability in real-world testing. 


II. Background and Literature Review 


2.1 Background Introduction 


This purpose of this chapter is to provide a basic background to the reader on 
geomagnetic mapping and modeling, geomagnetic measurements, and the filters pre- 
dominantly used in geomagnetic field navigation. Topics covered include the Interna- 
tional Geomagnetic Reference Field (IGRF) model of the Earth’s core field, external 
sources of the Earth’s magnetic field, and the magnetic anomaly field caused by the 
geology of the Earth’s crust. Scalar, vector and gradient measurements of the Earth's 
magnetic field will be introduced as well as the instruments used to take these mea- 
surements. The Kalman filter (KF) and extended Kalman filter (EKF) algorithms 


will be presented. Finally, a method for evaluating filter performance is described. 


2.2 Geomagnetic Mapping and Modeling 


The Earth's magnetic field, or the geomagnetic field, is comprised of the sum- 
mation of several individual magnetic fields. The source that accounts for about 
98% of the Earth's total geomagnetic field is the core field which is generated by 


electro-magnetic currents in the outer core of the Earth [12]. 


International Geomagnetic Reference Field. 


The IGRF is a spherical harmonic model of the Earth's core field. It is published 
every five years by the International Association of Geomagnetism and Aeronomy 
(IAGA) [31]. Measurements of the geomagnetic field collected from observatories, 
low-altitude aerial surveys, and satellite surveys provide the basic data for fitting 
the IGRF model [12]. This type of model approximately fits a periodic model onto 


a sphere with a set of coefficients [7]. The degree of the harmonics correspond to 


4 


a spatial wavelength [7]. Higher degree harmonics correspond to higher frequency 
signals and vice-versa. The IAGA is able to approximately isolate the magnetic field 
due to the core field of the Earth, because harmonics of the geomagnetic field of up 
to degree 13 (3,100 km wavelength and longer) are dominated by the core field [12] 
as shown in Figure 1. The IGRF models the core field by only including harmonics 
up to degree 13 [12]. The IGRF also includes approximations for secular variations, 
or long-term variations in both magnitude and direction of the core field, based on 


the actual rate of change from previous years [12]. 


F 


Core Field 


Log, (Power) 
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Crustal Field 
Temporal Variations 


0 8 16 24 32 40 48 56 64 72 
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Figure 1. Power Spectral Density of Earth’s Total Magnetic Field Spherical Harmonics 
[25] 


Temporal Variations. 


Variations in the total magnetic field that stem from sources external to the Earth 
are considered temporal variations. The strength of the magnetic field due to these 
external sources is weak compared to the magnetic field due to the Earth’s core. Also, 


temporal variations are often on a much shorter time scale than secular variations, 


so they can be approximately separated from the core field [25]. Because temporal 
variations occur in such rapid cycles, they are not captured by the IGRF. 

Temporal variations can be identified for navigation purposes by their frequency. 
Low-frequency variations tend to look like constant biases for short-term navigation, 
while the higher-frequency variations look like white noise [7]. The variations that 
occur at a similar frequency as data used for magnetic navigation are more difficult, 
if not impossible, to separate from the anomaly field [7]. 

The Earth’s ionosphere is a shell of ionized gas that reaches from about 50 km to 
beyond 1000 km above the Earth's surface. Movement of ions in the ionosphere result 
in electrical currents, which induce magnetic fields. One cause of the movement of 
these ions is solar heating [12]. Solar heating leads to current loops during the daytime 
hours, which causes a distortion in the Earth’s magnetic field as the Earth completes 
its daily rotation [31]. This cycle, as well as the gravitational pull of the moon, creates 
atmospheric tidal waves, which generally result in a variation of the total geomagnetic 
field of less than approximately 50 nanoTeslas (nT) over the course of a day, except 
during periods of increased solar wind [31]. 

Solar wind is radiation from the sun that causes a distortion in the Earth's core 
field into a comet-like shape that is known as the Earth’s magnetosphere as shown in 
Figure 2 [12]. Currents within the magnetosphere caused by the interaction of solar 
wind with the magnetic field are also a contributor to temporal variations. 

Solar storms are characterized as a period in time where the Earth is subjected to 
an unusually high amount of solar wind. The timing of solar storms is unpredictable 
and may lead to disturbances in the magnetic field of hundreds of nT or greater when 


the particles interact with the Earth’s magnetic field [12]. 
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Figure 2. Earth’s Magnetosphere [12] 


Crustal Sources. 


The final source of interest that contributes to the Earth’s total magnetic field is 
the Earth’s crust, which accounts for a small portion of the total field. The crustal 
field is made up of the superposition of all induced and remnant magnetization in the 
rocks and sediments that make up the Earth’s crust. Magnetization of a material in 
the Earth’s crust, a rock for example, consists of both induced and remnant compo- 
nents. The induced component of the magnetic field exists only in the presence of 
an external magnetic field, when the magnetic moment of the atoms that make up 
the material align under an external field, and adds to the total magnetization of the 
material [12]. Induced magnetization within a material is dependent upon the mag- 
netic susceptibility of that material. Magnetic susceptibility is the ease with which 
the material is magnetized by an external field [12]. 

Remnant magnetization is the component of a material’s total magnetization that 
is retained from a previous magnetic environment. The remnant component of the 


material’s magnetization is present without an external magnetic field and adds to the 


induced component to give the total magnetization of the material [12]. The direction 
of remnant magnetization is generally not in the same direction as the external field, 
as is the case with induced magnetization; However, the magnitude of a materials 
remnant magnetization versus its induced magnetization is generally much smaller 
[12]. 

The spatial variation of the crustal field varies at a high frequency relative to 
the core field. When flying over the crustal field, the frequency of the crustal field 
overlaps with the frequency of the temporal variations as described above, making it 


difficult to distinguish the two signals [7]. 


Magnetic Anomaly Mapping. 


Å magnetic anomaly is a vector deviation from a reference field. For the purposes 
of magnetic navigation, the reference field used is the Earth's core field, and the 
vector deviations, or magnetic anomaly vectors, come from the Earth's crustal field. 
The crustal field is relatively static and has a high spatial frequency, which makes 
it a good signal to use for magnetic navigation [7]. Maps have been created to 
capture a representation of the Earth's magnetic anomaly vectors, but a distinction 
must be made between the true anomaly vector at each observation point, and the 
scalar representation present in the magnetic anomaly maps at each corresponding 
observation point. 

The Earth’s total magnetic field vector By is approximately the vector sum of the 
Earth's core field (Bigrr) and the Earth’s crustal field (Ba) as shown in Equation 1. 
During magnetic surveys, care is taken to remove the effects of the temporal variations 
from the measurements, so the total field vector is an approximation of the sum of the 


Earth's core field and crustal field. Some of the effects from the temporal variations 


remain and are included in the total anomaly field vector [7]. 


B+ = Bicrr + Ba (1) 


Figure 3 shows a visual representation of this vector summation. B, is exaggerated 
in the figure for clarity and its intensity is more on the order of 2% of the intensity 


of Barr. 


BIGRF 


Figure 3. Earth's Total Field as a Vector Sum of Earth’s Core and Crustal Field 
Components 


Given current instrumentation limitations, only the intensity of the total field 
vector B; is usually collected during magnetic surveys as opposed to the intensity and 
direction of the vector. However, the intensity and direction of the core field vector 
(Bicre) is available from the IGRF model. Because we do not have the direction of 
the total field vector, we can not do a direct vector subtraction between the core field 
and the total field measurement to obtain the full anomaly vector shown in Figure 3. 
Instead, knowing that the intensity of the core field vector is much greater than the 
intensity of the anomaly field vector allows for the assumption that the anomaly field 
vector can not perturb the direction of the Earth's core field. Because of this, it can 
be assumed that subtracting the intensity of the core field vector from the intensity 
measurement of the total field vector gives å good approximation of the projection 
of the anomaly field vector in the direction of the core field. It is this approximate 
projection in the direction of the core field, or the amount that the anomaly vector 
"stretches" or "shrinks" the core field vector that is represented in the magnetic 


anomaly maps [7]. Figure 4 gives a visual representation of this assumption. 


Se 


— = > 
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Figure 4. Approximated Projection Presented in Magnetic Anomaly Maps 


Because |Bicrer| > |Ba], we assume the following: 


Buppros! = [Be] = |Bicrr| (2) 


2.3 Geomagnetic Measurements 


As described above, the Earth’s magnetic field is a vector field, so at any given 
point the direction of the magnetic field corresponds to the orientation of the vector, 


and the strength of the magnetic field corresponds to the length of the vector. 


Scalar Measurements. 


Scalar measurements of the magnetic field only capture the magnitude, or inten- 
sity, of the magnetic field. In Figure 5, the length of the vector B shown corresponds 
to the intensity of the magnetic field at that point. With an ideal scalar intensity 
sensor oriented in any direction at a given point, the scalar intensity measurement 
will be constant. 

Scalar intensity measurements include the magnitude of the Earth’s total magnetic 
field, so care must be taken to eliminate the effects of external fields and the core 


field to get the magnetic anomaly value. 


Vector Measurements. 


Vector measurements fully capture the intensity and direction of the magnetic 
field. The vector B shown in Figure 5 can be separated into its individual components 


B,, B, and B,. The magnitude of B is a scalar quantity and can be calculated with 
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Figure 5. Magnetic Field Vector Seperated Into Vector Components 


the following formula: 


B=/B2+B2+B! (3) 


When the length of Bz, By and B, are measured separately, we are able to calcu- 
late the intensity of the magnetic field vector as shown above as well as its direction. 
As with scalar intensity measurements, the vector B represents the Earth's total 
magnetic field including the magnetic field due to the Earth’s core, the Earth's crust, 
and external fields. The magnitude of each directional component of the magnetic 
field vector is different depending on the orientation of the sensor at any given point, 
even if the total intensity stays constant. When using a vector sensor as opposed to 
a scalar sensor, more care needs to be taken to accurately track the sensor’s attitude, 
because even the slightest error in attitude could lead to high errors when solving for 
the separate vector components. For example, if the green vector in Figure 6 is the 
true magnetic field vector, and we have 0.01 degrees of attitude error (resulting in the 
orange vector), we end up with 8.73 nT of error in the y-component of the magnetic 


field vector. 


11 


y B, = 50,000.00 nT 
By = 0.00 nT 
50,000 nT E: 
Í B, = 49,999.99 nT 
vs JE ee 
ME 
49,999.99 nT x 


Figure 6. Vector Components with and without Attitude Error 


With vector measurements, a coordinate frame transformation is required to trans- 
form the measurements from the body frame to the local or navigation (nav) frame. 
The body frame is aligned with the roll, pitch and yaw axes of the vehicle [39]. The 
x-axis points out the front of the vehicle, the y-axis points out the right side of the 
vehicle, and the z-axis points down. The nav frame is a local geographic frame that 
has its origin at the location of the vehicle. Its axes are aligned with the geographic 
north, east and down directions [39]. 

The magnetic field measurement from the sensor is in the body frame and is 


denoted as B, and is rotated to the nav frame using Equation 4 [27]. 


Ba = CBb (4) 


B, is the vector measurement expressed in the nav frame and CP is the direction 
cosine matrix required to transform the measurement from the body frame to the 
nav frame. The roll, pitch, and yaw angles (y, 6 and 4) of the vehicle may be used to 


calculate the transformation matrix, Cp. These calculations are shown in Equation 


12 


5 [13]. 


cos y cos Y + sin y sin W sind — cos ysin Y + sin y cos sind —sinycosd 
Cy = sin Y cos O cos W cos 0 sind 


sin y cos Y — cosysindsind —sinysinw—cosycosysin@ cosycosé 


(5) 


Gradient Tensor. 


Magnetic gradient measurements provide the difference between two magnetome- 
ter readings, either scalar or vector, taken simultaneously with a constant distance 
between them, or the spatial derivative of the magnetic field with respect to the 
distance between them. For our purposes, we will focuse on vector gradients. The 
partial derivative of the x-component of the magnetic field vector (B,) with respect 
to the x-direction is denoted as 0,B, or Bzz. If four sensors were configured on an 
axes as in Figure 7, Bar = By; — Bz3 , or the difference in the x-component of the 
magnetic field vector as the distance along the x-direction changes. Similarly, the 
partial derivative of the y-component of the magnetic field vector (B,) with respect 
to the x-direction is denoted as 0,B, or Bys. This would correspond to the difference 
between B,; and B,3 in Figure 7. Gathering these partial derivatives into a 3 x 3 


matrix produces the full nine-component magnetic gradient tensor, G. 


0, By Oy Ds -Bz Dos Bey Bzz 
aB 7 
= op = |%B, 0,B, 0,B,| = | Bye By By, (6) 
0, B; Oy By 0.B; Pye Bzy Bzz 
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Figure 7. Configuration of Four Vector Sensors 


where 


G is a symmetric and traceless matrix, meaning GT = G and Ba, + Buy + Bez =0 
[19]. Because the matrix is symmetric, we know that Bys = Bay , Ben = Baz and 
Bey = By, . Because the matrix is traceless, we know that Bz; = -(Bar + Byy) - 
These properties of the tensor leave us with five unique components out of the nine 
elements of the tensor: Baz, Bay, Bez, Byy, and Byz- 


As with vector measurements, the tensor measurements also require a transfor- 
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mation to put the local sensor frame measurements into the navigation frame. The 


equation for this transformation is below [27]: 


Ge OG 0 (7) 


where Gn is the tensor measurement rotated into the nav frame and Gy, is the body- 


frame measurement. 


Scalar Magnetometers. 


Different instruments are used to collect magnetic field data depending on which 
type of measurement is desired. For scalar measurements, nuclear resonance mag- 
netometers have sensors containing fluids or gases with properties that are sensitive 
to changes in the magnetic field. These are strictly scalar sensors that measure the 
absolute intensity of the magnetic field, but give no direction information [12]. The 
Geometrics airborne Cesium Resonance Magnetometer (G-823A) specifies a sensitiv- 
ity of 0.004 nT/vHz and absolute accuracy of less than 3 nT [11]. This high sensitivity 
and accuracy makes the Cesium Magnetometer a good example of an instrument that 


is used to collect absolute intensity measurements during magnetic surveys. 


Vector Magnetometers. 


Fluxgate magnetometers are one type of sensor used to measure the vector compo- 
nents of the relative magnetic fields [7]. The fluxgate within the sensor is a transducer 
that converts the magnetic field in one direction into a voltage. These sensors are 
small, durable, reliable and do not require much power to operate [23]. A Bartington 
Mag-03 Fluxgate magnetometer has similar sensitivity values to that of the Cesium 
Magnetometer (0.006 - 0.01 nT/vHz), but is much less accurate. With the Mag-03, 


if we are in an operating range of +70, 000 nT, we can get up to 350 nT of error given 


15 


the 0.5% scaling error listed in the data sheet [3]. Additional error is added when we 


take into account orthogonality errors in the mechanical sensor configuration. 


Magnetic Gradiometers. 


Using a vector gradiometer, as opposed to vector or scalar sensors as discussed 
above, virtually eliminates regional and temporal effects on the measurements [31]. 
This is due to the fact that the spatial gradient for temporal variations is nearly 
zero [7]. The Bartington Grad-13 consists of two Mag-03 Fluxgate Magnetometers 


separated by a distance of 1 m. This instrument is able to measure B,,, Byz, and Ba 


ya 
and may be visualized as sensor 1 and 3 from Figure 7. The Grad-13 has sensitivities 
similar to the scalar and vector sensors (0.02 nT /WHz/m), but with a higher accuracy 
than the vector sensor (+70 nT with the 0.1% scaling error listed in the data sheet) 
[2]. While the accuracy is far less than the scalar sensor, the gradiometer gives us 
much more information about the magnetic field at each location, to include Bz, By, 


and B, of the magnetic field at the reference sensor as well as Bzz, Byr; and Ba. 


yx 


Figure 8. Bartington Grad-13 Three-Axis Gradiometer [2] 


Using a Grad-13 gradiometer eliminates the need to perform the vector magne- 
tometer calibration as described above, because the Bartington gradiometer outputs 


corrected three-axis data from the vector-magnetometers [2]. 
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Å magnetic tensor configuration would allow us to record the same information as 
the magnetic gradiometer with the addition of the two unique gradient measurements 
required to calculate the full tensor matrix. The required configuration is similar to 
that displayed in Figure 7. Figure 9 shows five Honeywell HMR2300 Vector Magne- 
tometers arranged into a tensor configuration that was assembled at the Air Force 


Institute of Technology. 


Figure 9. AFIT Tensor Measurement Assembly 


2.4 Kalman Filter 


Kalman filtering is one method used to process noisy measurements such as magne- 
tometer measurements into optimal estimates of system random processes over time. 
In basic Kalman Filtering, the dynamics of a set of random variables can be mod- 


eled with Equation 8 for continuous-time processes or Equation 9 for discrete-time 
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processes. 


Xp+1 = OX, + Wr (9) 


The measurement function is shown in Equation 10. 


Zk = H,x;. + Vk (10) 


where 

x is the (n x 1) process state vector at time k 

$, is the (n x n) state transition matrix that relates the matrix x, to Xk+1 

B; is the (n x n) control input model matrix 

uz is the (n x 1) input matrix 

wx is the (n x 1) vector that includes the additive white Gaussian noise (AWGN) 
contribution to the state vector for the time interval (k + 1, k) 

F, is the (n x n) dynamics matrix 

zx is the (m x 1) measurement vector at time k 

H; is the (m x n) matrix that relates x; to Zk 

vi is the (m x 1) vector that includes the AWGN measurement error 


The covariance matrices for wą and v, are 


Elw,w;¿] = Qk (11) 


E[vavi ] = R; (12) 
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The error in the state estimate is defined as 


C= Xe Â, (13) 


where $, is the estimate of the state at time k. This estimate is based on all prior 
knowledge about the process prior to k. This is the estimate of the state before the 


measurement update has been applied at time k. 


P; = Ele; (e,)7] (14) 


The prior estimate X, is improved by implementing a measurement update as 
shown below: 


where x; is the updated estimate at time k and K, is the Kalman gain, which is the 
optimal gain that minimizes the mean-square estimation error and is computed by 


the following equation: 


K; = P; Hi (HP; HI +R,) (16) 


The error covariance matrix for the updated estimate is computed using the following 
formula: 


P, = (I- K,H,)P; (17) 


The updated estimate and the error covariance matrix are then propagated for- 
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ward in time using the formulas: 


Xp+1 = PX), (18) 


Pi = ®,Pı® + Qu (19) 


where Q, is the discretized covariance matrix for wą approximated by multiplying Q 
by the time difference between successive time steps. 

This process, which is considered the Kalman filter, continues recursively until 
estimates have been made for the entire trajectory. A complete derivation of the 


Kalman filter equations can be found in [6]. 


2.5 Extended Kalman Filter 


The Kalman filter algorithm is only optimal in systems that have linear dynamics 
and measurement models. When either the dynamics or measurement model are non- 
linear, the EKF may be used. The EKF assumes that the state transition function 
and measurement function are governed by nonlinear functions f and h respectively. 


The dynamics are modeled by the following equation: 
x, = f[x,, UL] + Wk (20) 
The measurement function follows the equation below: 
Zk = hixx| + Vk (21) 
Unlike in the Kalman filter, the matrix H} must be recalculated at each time 


step. Hy is the Jacobian of the measurement function h from equation 21 evaluated 
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at each time step. The measurement Jacobian is defined below: 


a Oh(x, k) 


Ox |x=; 


H, (22) 


The prior estimate at time k, X, is updated using equation 23 below, while the 


error covariance matrix is updated using Equation 24. 


P; = (I - K,¿H¿)P, (24) 
The Kalman gain is computed using Equation 25. 
K; = P; HI (HP, HT + Ry)! (25) 


Once the updated state estimate and error covariance matrix are calculated, they 


may be propagated forward using the following formulas: 


Xp = FX, un] (26) 


Pi = PP, Bi + Qa (27) 


where ®, is the matrix exponential of F,, which is the Jacobian of the function f 
from equation 20. Equation 28 defines the dynamics function Jacobian required to 


propagate the states and error covariance from time t; to time tr11- 


a Of(x, k) 


F 
i Ox x=%p 


(28) 


The state estimates and covariance are updated and propagated recursively until 


estimates have been made for the entire trajectory. Å detailed derivation of the EKF 
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algorithm may be found in [24]. 


2.6 Marginalized Particle Filter 


A specialized filter may be required to handle magnetic field navigation given the 
non-linearity of the map-matching problem. With inaccurate initial state estimates 
and no processing power limitations, the particle filter would be the filter of choice. 
The drawback of using the particle filter is its computational intensity. 

A particle filter is a simulation-based estimation technique used in filtering prob- 
lems to estimate the latent states of a complex dynamic statistical model, where after 
each observation, the state that gave rise to that observation is estimated [37]. 

The benefit of using a particle filter is its ability to handle highly non-linear 
dynamics and measurement models. As the number of particles approaches infinity, 
all possible values of each state can be estimated and the filter can determine the 
most likely state estimate. 

To estimate each state accurately, a large number of particles must be used. For 
a system model with nearly 20 states, the processing power required to perform 
the particle filter computations for all states is not viable [34]. To get around this 
limitation, an extension of the particle filter, the marginalized particle filter (MPF), 
may be used that separates the states that have linear dynamics from the states that 
have non-linear dynamics by partitioning the state vector. The linear states are then 
handled by a standard KF or EKF, while the non-linear states are handled by the 
standard particle filter. For a navigation filter, a majority of the states are linear, so 
the required processing power is reasonable when the MPF is employed. The MPF 


algorithm may be found in [34]. 
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2.7 Cramer-Rao Lower Bound 


The Cramer-Rao Lower Bound (CRLB) is the lowest possible state covariance of 
an estimator, such as the filter types described above [4]. [4] derived the CRLB for 
a terrain-aided navigation problem and concluded that the CRLB is equivalent to 
the error covariance in the EKF, with the measurement function h from Equation 
21 replaced with its gradient evaluated at the true state values at time k as opposed 
to the estimate of the state vector. This allows a comparison to be made between 
the CRLB and the actual filter error covariance to evaluate filter performance. Å 
high-performing filter should have error covariances close to, but not less than, the 
CRLB for all states. Additionally, the CRLB sets a lower limit on the DRMS of the 


filter and can be used to calculate how far from optimal the filter performs [4]. 


2.8 Background Conclusion 


In this chapter, an overview of the sources of the Earth's total magnetic field have 
been presented: the Earth's core field modeled in the IGRF, temporal variations due 
to external sources, and the Earth's crustal field. The magnetic anomaly field has 
been described as well as the distinction between the true anomaly field and the 
anomaly field value given in magnetic anomaly field maps. The difference between 
scalar, vector, and gradient measurements has been described and the instruments 
used to collect these measurements have been introduced. The Kalman filter and 
Extended Kalman filter algorithms were presented as well as a brief overview of the 
marginalized particle filter and an explanation of its benefits and limitations. The 


CRLB was presented to describe its merit in evaluating filter performance. 
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2.9 Literature Review 


Abundant research has been done on the field of non-GPS navigation, some of 
which pertains to magnetic field navigation. Navigation using the earth’s magnetic 
anomaly field has been investigated on several different platforms (aerial, sea-level, 
sub-sea-level, ground, and indoor). A majority of previous work for all platforms 
explored magnetic navigation using a scalar signal as opposed to the full magnetic 
field vector or tensor. However, other research has been done using the vector and 
tensor measurements of the magnetic field for UXO detection [17]. Each of these 
cases will be described in the paragraphs that follow. 

Research by A. Canciani [7] explored improving aerial navigation using the inten- 
sity of the magnetic anomaly field. Using this scalar measurement of the anomaly 
field, Canciani was able to demonstrate navigation accuracy of 13 meters distance 
root mean square (DRMS) in real flight tests. A limitation of Canciani’s research 
includes the use of one scalar measurement that is only able to capture the mag- 
nitude of the magnetic anomaly field. Additionally, high-accuracy results require a 
relatively high-velocity platform. The research may be expanded and improved upon 
using three measurements that capture all vector components of the anomaly field to 
increase navigation accuracy overall and at lower velocities (e.g. sea platforms). 

Research done by J. Wilson and R. Kline-Schoder [41] provides an aerial navi- 
gation solution using å magnetically-aided dead-reckoning system. This system uses 
air speed to provide å dead-reckoned navigation solution, with position updates pro- 
vided by the magnitude of the magnetic anomaly field. A tri-axial magnetometer 
was used, however, given the limitations of the magnetic anomaly map, these vector 
measurements were only used to calculate the scalar magnetic field strength. With 
this dead-reckoning system, position errors of 600 to 1200 m were achieved for flights 


after one hour with only moderate attempts to calibrate the magnetometers. Much 
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of this error could be due to aircraft noise, because the sensors were housed in the 
fuselage during flight. 

Three-axis magnetometer measurements and magnetic field maps were used in the 
dissertation by J. Shockley [36] to demonstrate self-contained ground-vehicle magnetic 
navigation. However, Shockley did not use inertial navigation system (INS) aiding 
to get high-precision results. Shockley was able to demonstrate navigation accuracy 
of 25-34 meters using a single magnetometer in a sedan and comparing it to a map 
he generated by sampling the magnetic field along the roadway at known positions. 
All measurements are taken in the vehicle's body frame and not converted to an 
absolute reference frame. Because of this, the mapping platform and the measure- 
ment platform are assumed to have the same attitude at each sample point and be 
traveling in the same direction between sample points. This sample point technique 
would not suffice for aerial applications where the measurements require a coordinate 
frame transformation. The results may be improved upon by pairing the three-axis 
magnetometer with an INS for higher-accuracy results that provide both a position 
and heading solution. 

W. Storms [38] used a three-axis magnetometer and a magnetic field intensity 
map for his research on indoor magnetic navigation. The heading angle was assumed 
to be known. The indoor navigation explored here required an INS and a previously 
generated map of the same indoor trajectory to be traversed by the user (taken at 
the same heading as the assumed user heading) for position updates. W. Storms was 
able to achieve position errors less than 0.2 m at all times with this method, but 
the solution only provided observability of the x- and y- position. This same method 
could be extended to estimate the heading as well. 

Similarly, research has been done in an indoor environment using four orthogonally 


mounted three-axis magnetometers to provide observability of the velocity of a person 
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traversing an indoor environment. The approach presented does not require the use 
of a magnetic map, because the trihedron gives a direct measurement of the magnetic 
field spatial partial derivatives [40]. 

D. Jeon fused odometer data with magnetic sensor data for indoor localization 
of an autonomous vehicle. An unscented Kalman filter (UKF) was used for the 
localization system that is based on magnetic markers in an indoor environment [14]. 
Å magnetic map was created by storing the location of magnetic markers. Å sensor 
was used to detect a magnet and match its reading to magnetic field data on the map 
to obtain an estimated position. The results show that the localization system was 
able to attenuate the cumulative position and heading errors of the vehicle that were 
present when only odometer measurements were used [14]. 

Sub-sea-level navigation presents its own unique challenges. Underwater platform 
navigation solutions frequently require INS aiding by surfacing to obtain a GPS fix. 
This is impractical and leaves vehicles such as military submarines vulnerable to de- 
tection. Previous research has been done by T. Karlsson [15] on terrain aided under- 
water navigation using a narrow-beam altimeter to aid the INS with a measurement 
of the depth directly below the vehicle. N. Kato [16] studied underwater naviga- 
tion of autonomous underwater vehicles using both geomagnetic measurements and 
bathymetric measurements. N. Kato concluded that underwater trajectories with a 
higher variation of bathymetric data gave higher accuracy navigation solutions. This 
is similar to aerial navigation solutions having higher accuracy when using a magnetic 
field map with a rich signal or traversing the magnetic field quickly leading to higher 
variation in magnetic field measurements [7]. Y. Huang presents a method to fully 
determine the attitude of an underwater vehicle to be used in underwater navigation 
problems [13]. The method is based on the use of full magnetic field tensor measure- 


ments. The algorithm worked successfully only when the initial angle error was less 
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than 20 deg and the sensor noise was less than 10 nT [13]. 

Minimal research has been done using magnetic underwater navigation given the 
challenges presented by the low velocity of å moving submarine, however the use of 
magnetic tensor gradiometers has been researched for underwater unexploded ord- 


nance (UXO) detection [17]. 
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III. Methodology 


Chapter III focuses on the detailed design of the navigation systems presented in 
this research. First, the reasoning for using the EKF in both navigation systems is 
presented. Then the detailed filter design for the navigation system that uses vector 
measurements is presented. This is followed by the filter design of the navigation 
system that uses tensor measurements. After the filter designs are laid out, the 


features of the simulation framework are outlined. 


Filter Type. 


The EKF was chosen for magnetic anomaly field navigation using vector mea- 
surements. In scalar magnetic navigation, the filter estimates a vehicle’s location by 
matching a magnetometer measurement to a location on the magnetic anomaly map 
that corresponds to the measurement. It is highly likely that there is more than one 
location with a map value closely matching the measurement. If a distribution of 
the possible map locations of the vehicle was generated after the filter received one 
measurement, this distribution of the latitude and longitude states would look highly 
non-Gaussian and would not fit a linearized model well. This is why the MPF was 
used for previous research into scalar magnetic anomaly field navigation [7]. 

Figure 10 shows the multi-variate distribution of the latitude and longitude states 
after a measurement update when a scalar measurement (B) was brought into the 
filter. When overlayed onto the magnetic field map, the shape of this multi-variate 
distribution is dependent on the shape of the contour line with a magnetic anomaly 
field value that matches the single measurement brought into the filter. 

The benefit of bringing three measurements into the filter as opposed to one is 
that the shape of the distribution is no longer conforming to the shape of a contour 


line, but of the intersection of three contour lines, one for each separate component of 
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the vector field. Figure 11 shows the distribution after a measurement update when 
three measurements (B,, By and B,) are brought into the filter. The distribution is 
overlayed onto the magnetic field maps of Bs, By and B}. This distribution looks 
more Gaussian than in Figure 10. An MPF is no longer required, because the EKF 
can, in theory, handle the Gaussian distribution that is created. This highly reduces 


the amount of processing power required. 
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Figure 10. Probability Distribution of the Horizontal Position States Overlaid on the 
Scalar Magnetic Anomaly Contour Map 


One drawback of the EKF in magnetic field navigation is its inability to handle 
multi-modal distributions. If the initial position uncertainty is large, the distribution 
of the latitude and longitude states after a measurement update could be multi-modal, 
leading to several possible position solutions. In order for the distribution to remain 
Gaussian and the EKF to perform well, it must be provided with accurate initial 
conditions. With an accurate estimate of the starting location, the distribution is not 
likely to be multi-modal. For this research, we operate under the assumption that the 
INS will be initialized using a true position from GPS, thus giving the filter accurate 


initial conditions. 
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Distribution of Horizonal Position States Given Vector Measurements 


Figure 11. Probability Distribution of the Horizontal Position States Overlaid on Com- 


ponent Magnetic Anomaly Contour Maps 


3.1 Navigation Filter Using Vector Measurements 


The first navigation system will use vector measurements of the anomaly field as 
well as inputs from an INS to aid the INS and constrain the position and angular 


drift within its sensors. In this section, the filter states will be introduced as well as 


the dynamics model and measurement function. 


Filter States. 


The filter estimates twenty states. The first seventeen states estimate INS errors 


and the last three states estimate the strength of temporal variations in each axis of 


the magnetic vector field. The state vector is following: 


Xk = [late lon, alt. v Uz v? Ex En Es hl 


€ 
T 
a. acct ac acc g? gY gi TV: TV» TVZ] 


where 
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lat., lon, and alt, are the INS position errors 


N „E 
e» Ves 


v and v? are the INS velocity errors in the north, east and down directions 
respectively 

€, €y and e, are the INS tilt errors about the x, y and z directions respectively 

h? is the aiding altitude error 

a. is the vertical acceleration error 

acc?, acc! and acc? are the INS accelerometer bias errors in the x, y and z axis 
respectively 

97, g? and gí are the INS gyroscope bias errors in the x, y and z axis respectively 


TV”, TV? and TV are the filter estimated temporal variations in the x, y and z 


components of the geomagnetic field respectively 


System Dynamics. 


The dynamics equation for the model used is a linearized function of the states. 
Because of the linearized dynamics in this system, the non-linear EKF dynamics 
function as shown in Equation 20 may be simplified to the linear dynamics function 
as defined in Equation 8. The dynamics matrix (F,.) represents a linearized 15-state 
Pinson Error model augmented with three states to model the temporal variations in 
each axis of the magnetic field vector, and two states to model the barometer errors. 

Variables that will appear in the dynamics matrix are following: 

lat is the INS solution for latitude in radians 

fn, fe, and fp are the north, east, and down specific forces from the INS 

Un, Ve, and vg are the north, east, and down velocities from the INS 

T, is the barometer error time constant 

kı, k2, and kz are barometer aiding constants used in the altitude aiding feedback 


loop 
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Ta is the accelerometer bias error time constant 

Tg is the gyroscope bias error time constant 

Try is the temporal variation error time constant 

re is the Earth's radius (6378135 m) 

w is the Earth's angular rate (7.2921151467 x 107 rad/sec) 

The dynamics matrix is shown below and the full derivation of the Pinson error 
model used may be found in [39]. 

ger (29) 
09. 11 F, 

Fp is the 11 x 11 Pinson model block. Fo is the block containing INS rotation 


matrices and F, represents the sensor and temporal variation bias error states. 
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The INS sensor errors and temporal variations are modeled as First Order Gauss 
Markov (FOGM) processes. A FOGM process X(t) has an exponential autocorrela- 


tion function that is defined in Equation 35|6]. 


Rx(r) = 02 All (35) 


where o and 7 are the standard deviation and time constant of the process respec- 
tively. 

The dynamics noise for this model is defined as 
2] = Q 


Elww 


Q = diag([0,,.3 VRW ıx3 ARW |x3 B 0 Aıxz G1x3 T1x31)20x20 (36) 


20? 

Bo EG 37 
2 (37) 
20? 

A, Se 38 
4 (38) 
20? 

G= = (39) 

g 

2 2 

ge (40) 
TTV 


where 
VRW is the noise strength of the velocity random walk in the INS 
ARW is the noise strength of the angular random walk in the INS 
op is the standard deviation of the barometer error 
dq is the standard deviation of the accelerometer error 
og is the standard deviation of the gyroscope error 


ory is the standard deviation of the temporal variations 
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B, A, G, and T are the driving noise strengths of the barometer, accelerometer, 
gyroscope, and temporal variation error states respectively 

For simplicity, the VRW, the ARW, the accelerometer, the gyroscope, and the 
temporal variations are assumed to have the same driving noise strength in each axis. 


This is why each of these appears three times in Q. 


Measurement Function. 


The measurement function is a non-linear function of the states as well as inputs 
from the INS. The measurement function used is the standard EKF measurement 
function as shown in Equation 41. Using an EKF allows for the linearization of the 
measurement function about the propagated estimate of each state’s trajectory. 


Equation 41 shows the measurement model for the body-frame measurements. 


Zk = h[xx] + Vk (41) 
Dex 
Li Byk (42) 
B, k 
4 Body Frame 


The measurement function, h, defines how the states relate to the magnetic anomaly 
field vector measurements and performs two main functions: map-matching and co- 
ordinate frame transformation. The coordinate-frame transformation is defined as a 
rotation from the nav frame to the body frame. 

The measurement function matches the INS reading of latitude and longitude 
(latrys and lonıns) combined with the states lat, and lon, at time tj. to the map of 
each vector component of the magnetic anomaly field. This is done using the mapping 
functions gz, gy, and g,. The mapping functions return the full vector measurement 


that we expect to see at the filter estimated location. We also expect to see the 
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effects of the temporal variations in the measurement, so the three TV, states are 
added to their respective vector components. This expected measurement is in the 
world frame, because the filter matches the location to nav-frame anomaly field vector 
maps. The nav-frame expected measurement must be rotated to the same frame as 
the actual vector sensor measurement. The expected measurement is rotated using 
the rotation matrix CP, which is the rotation matrix from the INS (CP ins) corrected 


using the tilt error states e,, €y, and €z. 


hrk Jallatıns, k + latek, lonıns, + loner) + TVE, 
hix] = Jaye] = Cex gy (latins, + later, ln in sx + loner) + TV}, (43) 
hz k ge(latıns, + latek, lonıns,r + loner) + TV, 


CP is the corrected rotation matrix calculated by applying a small angle correction 
to the rotation matrix from the INS (CP ins): The angle matrix e is used to solve for 
the rotation matrix, A, that satisfies Ch = CP ingA. 


The components of e at time k are the estimated tilt error angles, €x, €y, and €z. 


The magnitude of e is given by 


Eee (44) 


and 


lx =| e 0 eg (45) 


A is expressed in terms of the angle vector e in Equation 46 and a full derivation 


of this process may be found in [39]. 


ee O vik å 1 — cose 


Equation 47 applies the correction to the INS’ rotation matrix. 
CP = CP insA (47) 


The measurement Jacobian is defined as: 


Oh(x, k) 


ox |x=; 


H, ê (48) 

The partial derivative of each measurement with respect to every state must be 
calculated. The only states that show up in the measurement function are the hori- 
zontal position states (lateg, and lone p) within the mapping function, the tilt error 
states (€s, €y, and e,) within the small angle correction equation, and the tempo- 
ral variation states (TV%, TV%, and TV). The partial derivatives in the columns 


corresponding to the states that are not listed are zeros. 


Ohz,k Ohz,k Ohz,k 
latek lonek “°° OTV?, 
= Ohy k Ohy k Ohz,k 
H; latek lonek “°° OTV?, (49) 
Ohz k ðh- k ðh- k 
latek lonek ``! OTV?, 3x20 
= [Hato Hien, O3x4 He, He, He. 03x Hrv (50) 


3x20 


Equation 51 and 52 show how the first and second column of the measurement 
Jacobian were calculated. These columns contain the partial derivatives of h with 
respect to lat., and lon., respectively. There is no closed form solution for the 


spatial derivative of the maps within the f function, so they are obtained by finite 
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differencing. 


n,k 


where 


Gx (lat, + tat 
Jy (Lat, w Stat 

dlat 
gz(lat;, + “5, 


gx (lat, long 4 


gy (lat;,, long 4 


gz(laty, lon, 4 


lat, = latrns,k + late k 


lon; = lonıns,r + lone, 


long) — gz(latk — iat, long)) 
long) — gy(latk — Stat lon;)) (51) 
long) — gz (lat, — Y, long) 
- dl) — g,(laty, long — 32) 
- don) — gy(lat,, long — Hon) (52) 
- #2) — g.(latp, long — 4) 
(53) 
(54) 


The variables ólat and ólon can be chosen to be any small number. For the results 


presented in this research, a dlat and ólon of 0.5 m (converted to degrees latitude and 


degrees longitude respectively) is used. 


In Equation 50, the columns containing the partial derivatives of h with respect 


to the tilt states have a closed form solution. Each column was calculated separately, 


so will be presented below one column at a time. The column corresponding to the 
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€, State is shown in Equation 55. 


ge Je(latg, long) + TVS, 
H, = ghi = i (msa gy(laty, long) + TV}, (55) 
nee gz(laty,, long) TV}, 
gz(latgk, long) + T År 
= CP ing (SX) gy(lat,, long) + TV, (56) 
| g-(latg, long) + TVE, 


In order to differentiate the measurement function with respect to any of the tilt 


error states, the matrix CP in Equation 43 is replaced with CP ngA as shown in 


Equation 47, because A is a function of these states. In order to calculate oe from 
Equation 55, A was calculated symbolically in terms of the tilt error states as shown 
in Equations 44 through 46. A was then differentiated with respect to €s using a 
symbolic toolbox. 

The same process is followed to obtain the columns corresponding to the ey, and 


€z States. 


Oh, 3 
Oey 7 ge (lat, long) +7 sk 
í OAx 
_ | oh — cb 
H, = oe zu Cans ge, Jy (Lat, long) E TV); (97) 
oes gz(laty, long) + TV, 
er 92 (late, ons) + TV, 
OA 
_ | dh — cb k 
H. = ner = Cans ge, gy(latx, long) T TV; u 
Oh g-(latp, long) + TV}, 


(59) 


The columns of the measurement Jacobian corresponding to the three temporal vari- 
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ation states are shown in Equation 60. 


Ohz,k Ohz,k Ohz,k 
DTVå, OIV. OTVE, 


— | hye hyk ya | — b 
Hay OTVE, DVP, dek* Chk Loxa (60) 


Ohzk Ohz k Ohz,k 
OTV, DIVX, TVE, 


All of these columns are combined into the complete measurement Jacobian, Hj, 
as in Equation 50 to be evaluated each time step, k, given the current state estimates. 


Finally, the covariance of the white noise measurement error, v, is given by: 


Elvv*] = R3x3 (61) 


where the measurement noise strength for each measurement in z; is assumed to be 


equivalent. 


3.2 Navigation Filter Using Tensor Measurements 


The second navigation system will use tensor measurements of the anomaly field 
as well as inputs from an INS to estimate position and attitude of the vehicle. In 
this section, the filter states will be introduced as well as the dynamics model and 


measurement function. 


Filter States. 


Using tensor measurements allows for simplification of the dynamics model. The 
temporal variations were not included in the dynamics model because the gradient 
sensors used to measure the full tensor nearly eliminate any temporal and regional 
affects on the measurements. 


The tensor magnetic navigation filter estimates seventeen states as opposed to the 
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vector magnetic navigation filter’s twenty states. The states estimate INS errors and 
are equivalent to the first 17 states of the vector measurement state vector in section 


3.1. The state vector is 


Xz = [late lon, alte vN vi uP e, Ey €& hg 


€ € € € 


T 
di ac acc GE g? g? gi) 


System Dynamics. 


The same dynamics model is used in the tensor measurement filter with the ex- 
ception of the last three rows and columns of F in Equation 29, which correspond 
to the temporal variation states. These rows and columns are omitted in the tensor 
measurement dynamics matrix. 

As with the vector navigation model, the driving noise of the sensors are modeled 
as FOGM processes. The dynamics driving noise is the same as in Equation 36 with 
the exception of the last three rows and columns, which are omitted because the 


tensor measurement filter does not model the temporal variations. 


Measurement Function. 


The measurement function when using tensor measurements is of the same form 


as the EKF in Equation 21. The measurement vector zz is shown in Equation 62. 


Bree 
Byg k 

Zk = De (62) 
Byy,k 


Bzyk 
Body Frame 
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The measurement vector contains the five unique tensor components needed to 
completely describe the tensor. The filter requires the tensor measurements to be 
in full 3 x 3 matrix form at some points and 5 x 1 vector form at others. To deal 
with this transformation, the function M is defined in Equation 63. This function 


transforms a vector of size (5 x 1) to a symmetric, traceless matrix of size (3 x 3). 


ur = M (uy) (63) 


where 


ui 
U12 
Uy = |u3 (64) 
U22 
U23 
Uil Up U13 
UT = [uz Ua U23 (65) 


Us U —(u + U22) 


The inverse (M’) transforms a symmetric, traceless matrix of size (3 x 3) to a 


vector of size (5 x 1) containing its five unique components as shown in Equation 66. 


uy = M'(ur) (66) 
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where 


Uil U12 Us 
UT = [ua U22 U23 (67) 


U31 U32 U33 


Uy = |ui (68) 


As with the vector measurement function, the tensor measurement function per- 
forms the map-matching and coordinate frame transformation required. The measure- 
ment function matches each of the five unique tensor measurements to their respective 
maps using the mapping functions 922, Jys, Jzx Jyy, and gzy. These five mapping func- 
tions return the full tensor measurement that we expect to see at the filter estimated 
location in the nav frame. Unlike the vector measurements, the effect of the tempo- 
ral variations in the tensor measurements is negligible because the majority of their 
effects are canceled out when the spatial gradient is measured. Thus, we do not need 
to add the estimated effects of temporal variations into our expected measurement. 
The expected measurement returned from the mapping function is in the nav frame. 
It must be rotated into the same frame as the actual tensor sensor measurement as 
shown in Equation 7. The expected measurement is rotated using the rotation matrix 


CP which is calculated by using Equations 44 through 47. Equation 73 shows the full 
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tensor measurement function. 


lat, = latıns. + later (69) 
lon; = lonıns.r + lonek (70) 
9xx(laty, long) 
Gyx(latg, long) 
Sv (latr, long) = | g.. (lat;, long) (71) 
Jyy (Lat, long) 
Jzy (lat, long) 
Sr(lat;,, long) = M (gv (lat), lon;)) (72) 
Pez, k 
Ryx,k 
hie] = hr = | hank | = M (OR x gr (late, long) Cry.) (73) 
Nyy,k 
hzy,k 


The only states that show up in the measurement function used in this case are 
the horizontal position states (latex and lone g) within the mapping functions and the 


tilt error states (€,, €y, and e,) within the small angle correction equation. The partial 
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derivatives in the columns corresponding to states that are not listed are zeros. 


Ohgz,k hrr k hrr k 
Olatex  Oloner ``? Og? y 
Ohyak hyr k Ohyz,k 
Olat Olon nr Og? 

H, = ek ek Je k (74) 
Əhzy k Əhzy k Əhzy k 
latek lonek ``! Og k Bx17 

= Hat. Hion O3x4 He, He, He, 03x8 (75) 


5x17 


Hat. and Hion, are calculated using finite differencing as with the vector measure- 
ment Jacobian shown in Equations 51 and 52. 


H., H 


€x > €y > 


and H., are calculated using the product rule, which is defined in 
Equation 76, where q and w are arbitrary differentiable functions. Equations 78 


through 80 show the substitutions needed to find H., using the product rule. 


O(qw) ôq Ow 


de, J des | Tee, (76) 
m= = = ar (A) (77) 
qu = M(h,) = CP gr(lat,, long) ce.” (78) 
q= Chk = ChinsAk (79) 

w = gr (lat, long) Cho = gr(lat,, long) (Chis Ar)” (80) 


The functions q and w are then differentiated with respect to e, in Equations 81 and 


82 respectively and substituted back into Equation 76 to get the final H., in Equation 
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84. 


Og b OA; 

os MR G ER 81 
de n,INS dc. ( ) 
Ow OA; T 

de = gr(lat;, long) (Coms g) (82) 

O(qw OA OA 
(aw) = CP inser (late, long) (CP msAr)" + CP insAngr (late, long) (ns 
OE de, de, 

(83) 
(Uaw) 
BEE ( de ) (84) 
H., and H., are calculated using the same process: 

O(qw OA OA 
EE = CP ins gr(laty, long) (CP nsAr)" + CP insAngr (late, long) (msn 
de, Oey de, 

(85) 
,(OUqw 
HS ( s ») 20) 
y 

O(qw OA OA 
(aw) = C? ins Br (late, long) (Ch ns Ar)” =. CP insAnSr (lath, long) (ns 
De, De, Oe, 

(87) 
(qu 
By ae ( < ») (88) 
GAL, T and oat are calculated with a symbolic toolbox. 


All of these columns are combined into the the complete measurement Jacobian, 
H,, as in Equation 75, to be evaluated at each time step, k, given the current state 
estimates. 


The covariance of the white noise measurement error, v, is following: 


E[vv"] = Rsxs (89) 
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y 


y 


y 


where the measurement noise strength for each measurement in Zz; is assumed to be 


equivalent. 


3.3 Truth Trajectories 


Four different truth trajectories were used in the simulations. The different tra- 
jectories are used to demonstrate the difference in filter performance with varying 
vehicle velocities and positions. The actual boat trajectory used is characterized by 
low velocities and the actual aerial trajectory is characterized by high velocities. 

Additionally, the geographic location of the trajectory may affect filter perfor- 
mance. For example, the frequency content of the magnetic anomaly field at sea-level 
over deep seas is generally much lower than over shallow seas (e.g. over the continental 
shelf off the western coast of the U.S.) and may lead to degraded filter performance. 
To investigate the effect of geographic location on the navigation filter, the latitude 
and longitude of the truth trajectories may be shifted to alternate geographic loca- 
tions. This is possible because actual vector and tensor magnetic field measurements 
were not recorded during the course of the trajectories, so the flight and ship path 
are not tied to the actual geographic location. 


The trajectories included in the simulation are following: 


1. Actual 1 hour boat trajectory directly off the western coast of the U.S. 
2. Actual 1 hour aerial trajectory directly off the western coast of the U.S. 
3. Actual 1 hour aerial trajectory over the continental United States 


4. Actual 25 hour boat trajectory over the deep sea 


Figure 12 shows the geographic location of the four trajectories. Trajectory 1 and 
Trajectory 4 are both mapped by the orange path as Trajectory 1 is equivalent to the 


first hour of Trajectory 4. 
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Figure 12. Map of Trajectories 


The GPS position solution paired with the on-board sensors was used to generate 
the truth trajectories. While the truth data is not directly available to the navigation 
filter, it is used to calculate the amount of error present in the navigation filter solution 
and analyze overall filter performance. It is also used to generate simulated magnetic 


vector and tensor measurements as well as simulated INS solutions. 


3.4 Generating Vector and Tensor Maps 


The maps used for a majority of this research are segments from the North Amer- 
ican Magnetic Anomaly Database compiled by the U.S. Geological Survey (USGS). 
This is a database of aerial magnetic data collected and pieced together over the con- 
tinental United States and extending slightly over both coasts to cover shallow seas. 
[1]. Figure 13 shows the scalar anomaly field content collected over the entire North 


American continent. 


Figure 13. Scalar Magnetic Anomaly Content from NAMAD Data [1] 


The following figure shows a contour map of the scalar magnetic anomaly map 


used for the trajectories directly off the western coast of the United States. 
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Figure 14. Scalar Magnetic Intensity Contour Map from NAMAD Data 
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A Fourier method was used as in [5] to transform the scalar anomaly field into 
the vector components of the anomaly field. The individual directional components 


were found using the following equations: 


F [Bz] = HATE (a) (90) 
F[B,) = FIAT] F [yy] (91) 
F[B.] = FIAT] Fiy] (92) 
where 
ik; 
s 4 Ñ k ie pg ilke fr Tr Ku) sa 
ik 
Fipa) = — - 94 
Í k| fz + (Kafe ag ky fy) ) 
E IR 
= u u k b gg A a Rida) 


AT = total field anomaly measured in an ambient field 


(fa, A A) = unit vector in the direction of the ambient field 


fr = cos(I) cos(D) (96) 
fy = cos(I) sin( D) (97) 
fz =sin(1) (98) 


where F denotes the Fourier transform, k, k,, and ky are the wave numbers and I 
and D are the inclination and declination angles from the IGRF. 

The inverse Fourier Transform of Equations 90 through 92 gives the individual 
vector component maps used in the vector measurement navigation filter. These maps 


contain magnetic data in nT. 
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Figure 15. Contour Maps of Magnetic Anomaly Field Vector Components 


Calculating the spatial gradient of these vector component maps gives the five 


unique gradient maps of the full magnetic field tensor used in the tensor measure- 


ment navigation filter. These maps contain individual gradient values in nT/km and 


examples are shown in Figure 16. 
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Figure 16. Contour Maps of Unique Magnetic Anomaly Field Tensor Components 


Two simulation cases were run using maps generated from a source other than 


the NAMAD. For these cases the Enhanced Magnetic Model (EMM) was used to 
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evaluate the vector maps. Similar to the IGRF, the EMM is a spherical harmonic 
model. However, the EMM captures the magnetic anomaly field up to degree 790 as 
opposed to the IGRF capturing up to degree 13 as described in Chapter II [29]. The 
EMM is compiled from satellite, marine, aeromagnetic and ground magnetic surveys 
and resolves magnetic anomalies down to a 51 km wavelength [29]. Only capturing 
wavelengths down to 51 km leaves the EMM with a much lower resolution than the 
NAMAD. The EMM maps are not able to capture the high frequency content of the 
crustal field, which is what is primarily captured in the NAMAD maps. 

While the maps generated from the EMM are of a lower resolution, they model 
the magnetic field globally and are not limited to North America as the NAMAD 
data is. The EMM maps will be referred to as global maps for this reason. Figure 
17 shows the global EMM magnetic field data for the z directional component of the 
magnetic vector field. Vector component maps are directly resolved from the EMM. 
Figure 18 shows the contour maps of each vector field component generated from the 


EMM that were used for the simulation. 
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Figure 18. Contour Maps of Magnetic Anomaly Field Vector Components Resolved 


From the EMM 


The spatial gradients of these vector maps give maps of the five unique tensor 


components used in the simulations. Figure 19 shows contour maps of the five unique 


tensor components resolved from the EMM data that was used for the simulation. 
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NGDC-720 Version 3.0: Bz at Earth Surface 
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Figure 17. Vector Component Magnetic Field Content from EMM Data [29] 
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Figure 19. Contour Maps of Unique Magnetic Anomaly Field Tensor Components 


3.5 Generating Simulated Measurements 


Vector and Tensor measurements were not collected for the actual ship and aerial 
trajectories used. For the simulation, these measurements were generated using the 
process outlined below: 

1. To generate both vector and tensor measurements, uncorrupted magnetic field 

values are calculated by entering the truth trajectory values for latitude and 
longitude into the same mapping functions found in the respective measurement 


models. 


2. Vector measurements are corrupted with measurement biases that are modeled 
as FOGM processes to represent the magnetometer measurement error. Å o 
value of 3 nT and a 7 value of 600 seconds were chosen for this simulation 


framework. If tensor measurements are being generated, measurement biases 


are not added. 


3. These corrupted vector or tensor measurements are rotated into the body frame 


using the rotation matrix given in the truth trajectory. 


4. For both vector and tensor measurements, white Gaussian noise with a covari- 


ance of R = 0? is added. 


3.6 Generating Simulated INS Errors 


This simulation framework requires the INS solution as opposed to raw INS data. 
To create the realistic INS solution, a 17-state Pinson error model was used to gener- 
ate realistic INS errors throughout the entire truth trajectory. The truth trajectory 
position was then corrupted with the simulated INS position errors to get realistic INS 
position solutions. The true rotation matrix was also corrupted with the simulated 
INS tilt errors to get realistic INS solutions for rotation matrices. These simulated 
INS position and rotation matrix solutions are then used in the measurement func- 
tions of both navigation filters. Simulation parameters used to generate INS errors for 
a tactical-grade, navigation-grade, and 10X-grade INS are shown in Tables 1 through 
3. The 10X-grade INS model is characterized by having position and angular drift 
noise strengths (VRW and ARW) ten times less than the navigation-grade INS model 


noise strengths. 
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Table 1. Simulation Parameters Chosen for Tactical-Grade INS Cases 


Simulation Parameter | Chosen Value 

V RW, 0.3m/s?VHz 

ARW, 6.06 x 107 rad/sv Hz 
Og 7.27 x 1076 rad/s 

e 3600 s 

Ta 0.5 x 107° m/s? 

Ta 3600 s 

Tb 10m 

Tb 3600 s 


Table 2. Simulation Parameters Chosen for Navigation-Grade INS Cases 


Simulation Parameter | Chosen Value 
VRW, 0.001 m/s?V Hz 
ARW, 9.70 x 107? rad/sv Hz 
Og 1.45 x 1078 rad/s 

Ta 3600 s 

Oa 25 x 10 m/s? 

Ta 3600 s 

Tb 10m 

Ty 3600 s 
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Table 3. Simulation Parameters Chosen for 10X-Grade INS Cases 


Simulation Parameter | Chosen Value 

V RW, 0.0001 m/s?VHz 
ARW, 9.70 x 10° rad/sv Hz 
Tg 1.45 x 107°? rad/s 

Ti 3600 s 

Oa 25 x 1077 m/s? 

Ta 3600 s 

OD 10m 

Tb 3600 s 


Examples of the simulated position and angular drift for the tactical-grade INS 
modeled are shown below in Figures 20 through 23 respectively. These figures show 
the predicted standard deviation for the position and angular drift as well as examples 


of the error generated for use in the simulations. 
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Figure 20. Generated North INS Error Using Tactical-Grade INS Model 
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Figure 21. Generated Fast INS Error Using Tactical-Grade INS Model 
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Figure 22. Generated North Tilt INS Error Using Tactical-Grade INS Model 
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Figure 23. Generated Fast Tilt INS Error Using Tactical-Grade INS Model 


When the simulation parameters are changed to model the navigation-grade INS, 
the simulated position and angular drifts are less than when the tactical-grade INS 
model was used. The predicted standard deviation of these position and angular 
errors are shown in Figures 24 through 27. It is clear that when the navigation- 
grade INS was simulated, the estimated standard deviation of position drift grew to 
approximately 4.5 km within the hour, while the tactical-grade INS simulated position 
drift standard deviation reached almost 200 km. This pattern is consistent with the 
position and angular drift specifications of an actual tactical and navigation-grade 


INS [10]. 
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Figure 24. Generated North INS Error Using Navigation-Grade INS Model 


Predicted INS East Error Standard Deviation 


— East INS Error 
=== Predicted 1-0 


East Error (m) 


Time (min) 


Figure 25. Generated Fast INS Error Using Navigation-Grade INS Model 
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Figure 26. Generated North Tilt INS Error Using Navigation-Grade INS Model 
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Figure 27. Generated East Tilt INS Error Using Navigation-Grade INS Model 


63 


IV. Results 


The performance of the navigation systems introduced in the previous chapter are 
presented in Chapter IV. All simulation cases are listed and the metric used to evalu- 
ate filter performance is introduced. Each trajectory is then analyzed separately and 
the navigation accuracies and filter behavior for each navigation system combination 
are compared. Conclusions are drawn regarding the navigation system combination 


that gives the best navigation accuracies. 


4.1 Simulation Cases 


The specific cases investigated during this simulation case study are outlined in 
Table 4. Each case is characterized by the type of measurement, the truth trajectory, 
and the grade of INS used for navigation. Table 4 shows the specifics for each case. For 
simplicity, each specific case will be referred to by an abbreviation of the measurement 
type used within the navigation system (e.g. “VEC” for vector measurements and 
“TEN” for tensor measurements) combined with an abbreviation of the quality of 
INS used in the navigation system (e.g. “TACT” for tactical-grade INS and “NAV” 
for navigation-grade INS and *10X” for 10X-grade INS). For example, case 2 in Table 
4 below will be referred to as the VEC-NAV case for the coastal boat trajectory and 
case 7 will be referred to as the TEN-TACT case for the coastal airplane trajectory. 

For all cases, the initial position uncertainty was set to 200 m, the initial velocity 
uncertainty to 0.01 m/s and initial attitude uncertainty to 0.02 degrees. For all 
cases, a measurement noise covariance of R was used to generate measurements and 
a dynamics noise covariance of Q was used to generate INS data. For filter tuning, in 
the vector measurement cases, the filter’s measurement model was given a value of 2R 


and for the tensor measurement cases, the measurement model was given a value of 
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Table 4. List of Simulation Cases 


Case | Trajectory Measurement Type | INS Quality 
1 Coastal Boat Vector Tactical 

2 Coastal Boat Vector Navigation 
3 Coastal Boat Tensor Tactical 

4 Coastal Boat Tensor Navigation 
5 Coastal Aerial Vector Tactical 

6 Coastal Aerial Vector Navigation 
7 Coastal Aerial Tensor Tactical 

8 Coastal Aerial Tensor Navigation 
9 Continental Aerial | Vector Tactical 
10 Continental Aerial | Vector Navigation 
11 Continental Aerial | Tensor Tactical 
12 Continental Aerial | Tensor Navigation 
13 Deep-Sea Boat Vector 10X 

14 Deep-Sea Boat Tensor 10X 

15 Global Model Boat | Vector 10X 

16 Global Model Boat | Tensor 10X 


10R. In the vector measurement cases, the filter’s dynamics model was given a value 
of 2Q, while the tensor dynamics model was given a value of 10Q. This increased the 
stability of the filters, because the filter was expecting a noisier measurement than it 
was actually receiving. Setting the measurement covariances within the filter two or 
ten times greater than the covariance of the simulated measurements kept the filter 
from becoming over-confident in some cases. The tensor measurement filter required 
a greater multiple because bias error states were not modeled and the filter tended 


to become over-confident more often than in the vector measurement filter. 


4.2 DRMS Error 


The Distance Root Mean Square (DRMS) error is the metric used to measure 
filter performance for each simulation case. Monte Carlo simulations were run to 
evaluate the average DRMS error of the navigation filter in each case. This was used 


to compare filter performance between the different scenarios. The CRLB DRMS 
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Error was also calculated to get an idea of how the navigation filter is performing 
compared to its own theoretical optimal performance. Equations 99 and 100 define 


the EKF DRMS error and CRLB DRMS error calculations respectively [28]. 
> her (de)? 


DRMS Error = (| S=—— (99) 
n 


On + Op) 
n 


CRLB DRMS Error = y Dial (100) 


where 
n = number of time steps in the trajectory 
dą = Euclidian distance between horizontal truth position and filter solution 


on and og are the CRLB standard deviations for the horizontal position states 


4.3 Coastal Boat Trajectory Results 


This section outlines the simulation results for the cases 1 through 4 that used 
the boat’s truth trajectory off the west coast of the U.S. The trajectory was one hour 
long and had an average velocity of approximately 6 m/s. 

Figure 28 shows the path of the ship over the three vector component maps of the 


magnetic field. 
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Figure 28. Boat Trajectory Overlayed Onto Vector Field Maps 


An example of the corrupted vector and tensor measurements generated for this 


trajectory are shown in Figures 29 and 30 respectively. The measurements were 


corrupted with zero-mean White Gaussian noise with a covariance of 1 nT? for the 


2 
vector measurements and 0.1 E for the tensor measurements. 


Table 5. DRMS Results for U.S. Western Coast Boat Trajectory 


Case | Measurement Type | INS Quality | DRMS Error | CRLB DRMS Error 
1 Vector Tact 420.54 m 342.39 m 

2 Vector Nav 185.00 m 241.39 m 

3 Tensor Tact 755.48 m 136.21 m 

4 Tensor Nav 35.94 m 49.36 m 


Table 6. Filter Error in North and East Tilt Error States for U.S. Western Coast Boat 


Trajectory 
Case North Tilt Filter Error | East Tilt Filter Error 
VEC-TACT | 0.07 deg 0.11 deg 
VEC-NAV |9x10°? deg 3 x1073 deg 
TEN-TACT | Unstable Unstable 
TEN-NAV | 1.2 x107 deg 1.9 x10 deg 
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Figure 29. Corrupted World-Frame Vector Magnetometer Measurements for the VEC- 
TACT and VEC-NAV Cases 
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Figure 30. Corrupted World-Frame Magnetic Tensor Measurements for the TEN- 
TACT and TEN-NAV Cases 
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CRLB Trend. 


The CRLB DRMS errors shown in Table 5 display the trend that should be 
expected from the EKF DRMS Error results. The VEC-NAV case out-performed 
the VEC-TACT case with 101 m less position error and the TEN-NAV case out- 
performed TEN-TACT with 86.85 m less position error. From this we can expect 
to see the navigation filter paired with a navigation-grade INS to out-perform the 
same filter paired with a tactical-grade INS. Additionally, the TEN-TACT and TEN- 
NAV cases, where tensor measurements were used, out-performed the VEC-TACT 


and VEC-NAV cases, where vector measurements were used. 


Monte Carlo Simulation Results. 


Figures 31 through 34 show the results of a 500 run Monte Carlo simulation for the 
VEC-TACT and VEC-NAV cases in the horizontal position states. Results from the 
Monte Carlo simulations for cases VEC-TACT and VEC-NAV reflect the expected 
trend from the CRLB results. With a DRMS Error of 185.00 m, the filter paired 
with the navigation-grade INS in the VEC-NAV case out-performed the filter paired 
with the tactical-grade INS in the VEC-TACT case by 235.54 m. Both filters were 
stable, with no divergent runs. While the position error in the VEC-TACT case 
increases to an average of 1.5 km approximately thirty minutes into the trajectory, 
the filter is able to lock back down on position and does not continue to drift as the 
un-aided INS would. In the VEC-NAV case, we do not see the same increase in error 
at 30 minutes, and the filter remains locked on to the true position with steady error 
throughout the entire trajectory. Figures 31 through 34 show the standard deviation 
of the Monte Carlo error, which tends to be slightly lower than the CRLB. The CRLB 
is the theoretical lowest possible state covariance, however, because of the tuning that 


was described previously (the filters being given a multiple of R and Q), it is possible 
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for Monte Carlo DRMS error to be slightly lower than the DRMS error of the CRLB. 


The CRLB is expecting greater measurement noise than is actually present in the 


corrupted measurements. The CRLB and the filter models are both given the same 


tuning parameters, so while the filter was able to out-perform the CRLB DRMS error 


in simulation, the filter predicted covariance matches, but does not exceed, the CRLB. 
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Figure 31. EKF Error in North Position State for the VEC-TACT Case 
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CRLB to EKF Comparison for East Filter Error 
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Figure 32. EKF Error in East Position State for the VEC-TACT Case 
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Figure 33. EKF Error in North Position State for the VEC-NAV Case 
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CRLB to ld Comparison for Fast Filter Error 


150 pri CMR ONT GT 
7 ee 
| ay 


ih 
MA 
PGA 


Fu NN 
IN a TN iy 


East Error (m) 


0 10 20 30 40 50 60 
Time (min) 


Figure 34. EKF Error in East Position State for the VEC-NAV Case 


The filter predicted covariance more closely matches the CRLB, as shown in Figure 
35 and 36. The DRMS error for the single runs plotted in Figure 35 and 36 were 
238.25 m and 152.64 m respectively. 
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Figure 35. EKF Error in Horizontal Position States for a Single Run - VEC-TACT 
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North Filter Error for a Single Filter Run 
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Figure 36. EKF Error in Horizontal Position States for a Single Run - VEC-NAV 


While the TEN-TACT and TEN-NAV cases reflect this same trend of the navigation- 
grade INS out-performing the tactical grade INS, the TEN-TACT case, where tensor 
measurements were used demonstrated a DRMS error of 755.48 m. This error is 
greater than the filter in the VEC-TACT case, which also uses a tactical grade INS, 
but was only using three vector measurements as opposed to five tensor measure- 
ments. This does not agree with the CRLB DRMS error results. While the filter in 
the TEN-TACT case did not diverge, the position error plots in Figures 37 and 38 
display abnormal filter behavior and may be due to filter tuning. 

This same trend occurred for the TEN-TACT case during the coastal and con- 
tinental aerial trajectories as well. The tensor measurement navigation filter often 
diverged in Monte Carlo simulation when paired with a tactical grade INS. This was 
not a problem when the navigation-grade INS was paired with the tensor measure- 
ment filter for all TEN-NAV cases. The Monte Carlo position error plots for the 
TEN-TACT and TEN-NAV cases during the coastal boat trajectory are shown in 
Figures 37 through 40. 
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Figure 37. EKF Error in North Position State for TEN-TACT Case 
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Figure 38. EKF Error in East Position State for TEN-TACT Case 


74 


ÅR 


- Filter Runs 


Filter MC Mean +/- STD 


North Error (m) 


Time (min) 


Figure 39. EKF Error in North Position State for TEN-NAV Case 
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Figure 40. EKF Error in East Position State for TEN-NAV Case 


The navigation filter in the VEC-NAV and TEN-NAV cases were able to greatly 
improve upon the unaided navigation-grade INS position error drift that we would 
expect to be approximately 1 nmi/h (or 1,852 meters per hour) [10]. Additionally, 
the VEC-TACT case was able to improve drastically upon the drift expected from an 
unaided tactical-grade INS (approximately 10 nmi/h [10]). 

Aside from the accuracy of the position solution, the ability of the filter to accu- 
rately estimate the tilt error states is also examined. Their use in the measurement 
function to correct the INS rotation matrix makes their accuracy imperative. Accu- 
rate estimates of the tilt error states allow for higher accuracy resolution of the roll, 
pitch, and yaw of the vehicle throughout the trajectory. 

The tilt error states for the TEN-NAV case are shown for a single filter run in 
Figures 41 through 43. The DRMS error for this single run was 40.15 m. 
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Figure 41. EKF Error in North Tilt Error States for the TEN-NAV Case 


76 


East Tilt Filter Error for a Single Filter Run 
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Figure 42. EKF Error in East Tilt Error States for the TEN-NAV Case 
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Figure 43. EKF Error in Down Tilt Error States for the TEN-NAV Case 
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For the VEC-NAV case, where vector measurements were used, the tilt error states 
displayed the same stability as in the TEN-NAV case. The steady state values of filter 
error for the north and east tilt error states are listed in Table 6. 

When a tactical grade INS was used, it is clear that the tilt error state estimates are 
stable, but are less accurate overall than in the VEC-NAV case. The tilt error states 
for the single run of the VEC-TACT case are shown in Figure 44. The tactical-grade 
INS did not give as much observability of the tilt error states as the navigation-grade 


INS. 
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Figure 44. EKF Error in Tilt Error States for the VEC-TACT Case 


4.4 Coastal Airplane Trajectory Results 


The trajectory used for cases 5 through 8 was a one-hour airplane trajectory with 
an average velocity of 62.85 m/s. Figure 45 shows the path of the airplane over the 
three vector component maps of the magnetic field. An example of the vector and 


tensor measurements generated for these cases are shown in Figures 46 and 47. 
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Figure 45. Airplane Trajectory Overlayed Onto Vector Field Maps 
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Figure 46. Corrupted World-Frame Vector Magnetometer Measurements for the VEC- 
TACT and VEC-NAV Cases 
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Figure 47. Corrupted World-Frame Magnetic Tensor Measurements for the TEN- 
TACT and TEN-NAV Cases 


Table 7. DRMS Results for U.S. West Coast Aerial Trajectory 


Case | Measurement Type | INS Quality DRMS Error CRLB DRMS Error 
5 Vector Tact 254.96 m 315.26 m 

6 Vector Nav 103.43 m 128.18 m 

7 Tensor Tact Filter Diverged | 252.46 m 

8 Tensor Nav 61.07 m 85.97 m 


Table 8. Filter Error in North and East Tilt Error States for U.S. West Coast Aerial 


Trajectory 
Case North Tilt Filter Error | East Tilt Filter Error 
VEC-TACT | 0.07 deg 0.1 deg 
VEC-NAV | 1 x10 deg 1 x10 deg 
TEN-TACT | Unstable Unstable 
TEN-NAV | 6 x10* deg 1 x107 deg 


The CRLB results show that we expect the TEN-TACT case to have improved 


filter performance over the VEC-TACT case. And similarly, we expect to see the 
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TEN-NAV case perform better than the VEC-NAV case. The CRLB DRMS Error 
for the filters paired with a navigation-grade INS is lower than the same filters paired 
with a tactical-grade INS. 

As mentioned previously, the tensor measurement filter was unstable when a 
tactical-grade INS was used and the EKF DRMS results for the TEN-TACT case 
reflect this. However, the vector measurement filter paired with the tactical-grade 
INS remained stable for a majority of the runs (0.4% of filter runs diverged). Figures 
48 and 49 show the filter error in the horizontal position states for a 500-run Monte 


Carlo simulation of the VEC-TACT case. 
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Figure 48. EKF Error in North Position State for the VEC-TACT Case 
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Figure 49. EKF Error in East Position State for the VEC-TACT Case 


With the instability of the filter in the TEN-TACT case, the mean error of the 
Monte Carlo runs appears to grow without bound, and is considered divergent. Ås 
with the VEC-NAV and TEN-NAV cases during the coastal boat trajectory described 
earlier, the VEC-NAV and TEN-NAV cases during the coastal aerial trajectory also 
show stable filter performance. The Monte Carlo simulation results for the VEC- 
NAV and TEN-NAV cases during the coastal aerial trajectory are shown in Figures 
50 through 53. 
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Figure 50. EKF Error in North Position State for the VEC-NAV Case 
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Figure 51. EKF Error in East Position State for the VEC-NAV Case 
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CRLB to EKF Comparison for North Filter Error 
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Figure 52. EKF Error in North Position State for the TEN-NAV Case 
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Figure 53. EKF Error in East Position State for the TEN-NAV Case 
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The tilt error states for the TEN-TACT case, where a majority of the Monte 
Carlo runs diverged, remained at the CRLB until approximately thirty minutes into 
the trajectory, where the filter error began to grow unbounded. 

Changing the grade of INS to navigation (as in the TEN-NAV case) as opposed 
to tactical, we were able to achieve filter stability. The Monte Carlo simulation error 


for the north and east tilt error states for the TEN-NAV case are shown below. 
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Figure 54. EKF Error in North Tilt Error States for the TEN-NAV Case 
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Figure 55. EKF Error in East Tilt Error States for the TEN-NAV Case 


4.5 Continental Airplane Trajectory Results 


The trajectory from cases 5 through 8 was moved over land as shown in Figure 12 
for use in cases 9 through 12. Figure 56 shows the path of the airplane over the three 
vector component maps of the magnetic field. Examples of measurements generated 
at this location are shown in Figures 57 and 58. While the vector measurements are å 
steadily varying signal over the course of the trajectory, the tensor components show 
extreme variations for the first ten minutes of the trajectory and settle out for the 
remainder. Ås the airplane moved northwest as shown in Figure 56, it is clear where 


the high frequency content of the magnetic field drops down to low frequencies. 
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Figure 56. Airplane Trajectory Overlayed Onto Vector Field Maps 
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Figure 57. Corrupted World-Frame Vector Magnetometer Measurements for the VEC- 
TACT and VEC-NAV Cases 
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Figure 58. Corrupted World-Frame Magnetic Tensor Measurements for the TEN- 


TACT and TEN-NAV Cases 


Table 9. DRMS Results for U.S. Continental Aerial Trajectory 


Case | Measurement Type | INS Quality | DRMS Error | CRLB DRMS Error 
9 Vector Tact 514.53 m 549.13 m 

10 Vector Nav 141.17 m 182.7 m 

11 Tensor Tact Filter Diverged | 502.85 m 

12 Tensor Nav 205.95 m 136.24 m 


Table 10. Filter Error in North and East Tilt Error States for U.S. Continental Aerial 
Trajectory 


Case North Tilt Filter Error East Tilt Filter Error 
VEC-TACT | 0.12 deg 0.1 deg 

VEC-NAV |8x10°? deg 1.4 x10 3 deg 
TEN-TACT | Unstable Unstable 

TEN-NAV | 8x10°? deg 3 x10 deg 


The tensor measurement filter paired with a tactical-grade INS (the TEN-TACT 


case) remained as unstable as with previous trajectories. Looking closely at when the 
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filter diverges leads us to believe that the instability of the filter is compounded by 
the fact that the measurements remain fairly constant and unchanging after the ten 
minute mark. 

Figure 59 shows the filter error for the horizontal position states during a single 
filter run of the TEN-TACT case as well as the filter predicted covariance. The filter 
predicted covariance becomes low (over-confident) during the time that the tensor 
measurements are reflecting the rich signal at ten minutes. The filter confidently locks 
onto an incorrect position solution at this point judging by the steadily increasing 


error and continuing overly-confident covariance. 
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Figure 59. EKF Error in Horizontal Position States for a Single Run - TEN-TACT 


Case 


When a navigation-grade INS was used in the TEN-NAV case, the position error 
decreased and the filter was relatively stable during the Monte Carlo simulation. The 


Monte Carlo error for the horizontal position states for the TEN-NAV case are shown 


in Figures 60 and 61. 
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Figure 60. EKF Error in North Position State for the TEN-NAV Case 
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Figure 61. EKF Error in East Position State for the TEN-NAV Case 
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For a single run with a navigation-grade INS, the TEN-NAV filter had the same 
decrease in filter-predicted covariance right before ten minutes, but became less confi- 
dent as the measurements became less rich. The filter did not lock onto the incorrect 
position solution because of its over-confidence as it did in the TEN-TACT case. 

The behavior of the tilt error states in the TEN-TACT and TEN-NAV cases 
parallel the behavior of their position states. Figures 62 shows the filter error in the 
tilt error states and filter predicted covariance for a single run of the TEN-TACT 
case while Figures 63 through 65 show the filter error in the tilt error states for the 
TEN-NAV case. In the TEN-TACT case, the filter incorrectly locked on to a tilt 
error solution, whereas in the TEN-NAV case, the filter was able to clamp down on 
the correct north and east tilt error solution. For the down tilt error state, the filter 
estimate shifted while the measurements were rich and the filter was confident. The 


estimate stayed locked on to the same incorrect solution for the rest of the trajectory. 
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Figure 62. EKF Error in Tilt Error States for a Single Run - TEN-TACT Case 
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Figure 63. EKF Error in North Tilt Error States for a Single Run - TEN-NAV Case 
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Figure 64. EKF Error in East Tilt Error States for a Single Run - TEN-NAV Case 
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Down Tilt Filter Error for a Single Filter Run 
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Figure 65. EKF Error in Down Tilt Error States for a Single Run - TEN-NAV Case 


The VEC-TACT and VEC-NAV cases behaved similarly during the continental 
airplane trajectory as they did during the coastal airplane trajectory. The VEC- 
TACT filter during the continental airplane trajectory was relatively stable with 
only 4% of filter runs diverging, but was improved upon by pairing the filter with 
a navigation-grade INS. When the vector measurement filter was paired with the 
nav-grade INS in the VEC-NAV case, it was able to out-perform the tensor measure- 
ment filter the TEN-NAV case. The VEC-NAV case had a Monte Carlo simulation 
DRMS Error of 65.78 m less than the TEN-NAV case. This is most likely a result 
of the tensor measurements behavior around ten minutes. The Monte Carlo error for 
position states and tilt error states in the VEC-NAV case are shown in Figures 66 


through 68. 
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Figure 66. EKF Error in North Position State for the VEC-NAV Case 
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Figure 67. EKF Error in East Position State for the VEC-NAV Case 
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CRLB to EKF Comparison for North Tilt Filter Error 
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Figure 68. EKF Error in Tilt Error States for the VEC-NAV Case 


4.6 Deep-Sea Boat Trajectory Results 


The trajectory used in cases 13 and 14 was a 25 hour long boat trajectory in the 
same location as the coastal boat trajectory used in cases 1 through 4. The trajectory 
remained the same, but the maps used to navigate were altered to be more realistic 
for a location with an ocean-depth of 6000 m. A Fourier method was used to upward 
continue the anomaly map to an altitude of 6000 m. At this height, the frequency 
content of the map is much lower [5], and we expect to see a signal that is not as ideal 
for magnetic navigation. This high-altitude map was used to generate the vector and 
tensor maps via the Fourier method described previously. 

For both cases (VEC-10X and TEN-10X), a 10X-grade INS was used. A military 
ship navigating for longer trajectories (such as the 25 hour trajectory used here) would 
likely have a high-quality INS on board, and we model a higher-quality INS with the 
10X-grade INS as described in Table 3. This INS is modeled to have a position error 


of approximately 1 nmi after 24 hours [10]. 


Figure 69 shows the path of the ship over the three vector component maps of 
the magnetic field. An example of the corrupted vector and tensor measurements 
generated for this trajectory are shown in Figures 70 and 71 respectively. The mea- 
surements were corrupted with zero-mean White Gaussian noise with a covariance of 
27 nT? for the vector case and 0.8 u for the tensor case. The noise strengths were 
chosen to triple the measurement error of all previous cases. Greater measurement 
noise is expected in these cases, given the possible effects of the deep ocean on the 


Earth’s magnetic field at sea-level. 
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Figure 69. Boat Trajectory Overlayed Onto Vector Field Maps 
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Figure 70. Corrupted World-Frame Vector Magnetometer Measurements for the VEC- 
10X Case 


u Ensemble Corrupted Magnetometer Measurements (Zero Mean) 
T i T T 


40 


20 


j | | 


Ki m Y AM Ml EM VA 
ey Ay fol pana, Apr Na 
Vi | 


-60 i | | i 
0 


-20 


Magnetic Gradient (nT/km) 
o 


-40 


Time (hours) 


Figure 71. Corrupted World-Frame Magnetic Tensor Measurements for the TEN-10X 


Case 
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Table 11. DRMS Results for Deep-Sea Boat Trajectory 


Case | Measurement Type | INS Quality | DRMS Error | CRLB DRMS Error 
13 Vector 10X 194.90 m 249.23 m 
14 Tensor 10X 185.52 m 169.21 m 


Table 12. Filter Error in North and East Tilt Error States for Deep-Sea Boat Trajectory 


Case North Tilt Filter Error | East Tilt Filter Error 
VEC-10X | 2.6 x107% deg 3 x10* deg 
TEN-10X | 3.4 x10"* deg 4 x10* deg 


From Table 11, it is clear that both navigation systems were able to cut the drift of 


the INS down from the expected 1 nmi over the 25 hour trajectory to approximately 


200 m. The TEN-10X case only had approximately 10 m less DRMS error than the 


VEC-10X case for the 500 run Monte Carlo simulation. The position error states for 


a single run of the TEN-10X case are shown below: 
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Figure 72. EKF Error in North Position State for a Single Run - TEN-10X Case 
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Figure 73. EKF Error in East Position State for a Single Run - TEN-10X Case 


In both the VEC-10X and TEN-10X cases, the tilt error states remained minimal. 
Throughout the entire 25 hour trajectory, using these navigation systems, we would 
have an accurate attitude solution in addition to the relatively low DRMS errors. 
The error plots for the tilt error states in a single filter run of the TEN-10X case are 
shown in Figure 74 through 76. These plots show that the TEN-10X filter was able 


to retain an accurate orientation solution throughout the entire trajectory. 
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Figure 74. EKF Error in North Tilt Error State for the TEN-10X Case 
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Figure 75. EKF Error in East Tilt Error State for the TEN-10X Case 
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Down Tilt Filter Error for a Single Filter Run 
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Figure 76. EKF Error in Down Tilt Error State for the TEN-10X Case 


4.7 Boat Trajectory Results Using a Global Model 


The 25-hour boat trajectory was also used for cases 15 and 16. With the deep- 
sea boat trajectory above, the navigation systems were paired with simulated lower- 
frequency maps and were still able to perform. For the next two cases, the VEC-10X 
and TEN-10X navigation systems were paired with maps generated from the global 
EMM model. The trajectory used in these cases was shifted over an area of the 
continental United States to match the EMM map data available at the time of the 
simulation. These cases demonstrated the navigation accuracies possible for both the 
VEC-10X and TEN-10X navigation systems paired with low-frequency global model 
maps. 

Figure 77 shows the path of the ship over the three vector component maps of 
the magnetic field. Examples of corrupted vector and tensor measurements generated 


for the 25-hour boat trajectory using the global model maps are shown in Figures 78 


101 


and 79 respectively. The measurements for these two cases were also corrupted with 
white Gaussian noise with a covariance of 27 nT? for the vector case and 0.8 aT for 


the tensor case. 
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Figure 78. Corrupted World-Frame Vector Magnetometer Measurements for the VEC- 
10X Case 
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Figure 79. Corrupted World-Frame Magnetic Tensor Measurements for the TEN-10X 


Case 


Table 13. DRMS Results for Boat Trajectory Using a Global Model 


Case | Measurement Type | INS Quality | DRMS Error | CRLB DRMS Error 


15 Vector 10X 219.28 m 275.9 m 


16 Tensor 10X 38.10 m 94.98 m 


Table 14. Filter Error in North and East Tilt Error States for Boat Trajectory Using 
a Global Model 


Case North Tilt Filter Error | East Tilt Filter Error 
VEC-10X | 4 x10* deg 5 x10* deg 
TEN-10X | 8.8 x10”° deg 2 x10"* deg 


From Tables 13 and 14, it is clear that both navigation systems improved upon the 
expected unaided drift of the INS. The TEN-10X filter was able to achieve navigation 
accuracies of 38 m using the global map and steady state tilt errors down to 8.8 x 1073 
degrees. Figures 80 and 81 below show the north and east filter error for a single run 


using the TEN-10X navigation system. 
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Figure 80. EKF Error in North Position State for the TEN-10X Case 
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Figure 81. EKF Error in East Position State for the TEN-10X Case 


The error in the north and east tilt error states for this same filter run are shown 


below in Figures 82 and 83. 
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Figure 82. EKF Error in North Tilt Error State for the TEN-10X Case 
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Figure 83. EKF Error in East Tilt Error State for the TEN-10X Case 


These are promising results given the global availability of the EMM-720 model. 
These cases demonstrate that we are able to navigate anywhere around the globe 


using the navigation systems presented in this research paired with low-resolution 
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global model maps. When the VEC-10X and TEN-10X navigation systems are used 


in simulation, we are not limited by the low resolution of the available global model. 


4.8 DRMS Error for all Cases 


Table 15 lists the Monte Carlo simulation DRMS error for each case. It reflects 
the results presented above, but allows for comparison between the navigation system 


results for different trajectories. 


Table 15. DRMS Error Results for all Simulation Cases 


Case | Trajectory Meas Type | INS Quality | DRMS Error 
1 Coastal Boat Vector Tactical 420.54 m 

2 Coastal Boat Vector Navigation 185.00 m 

3 Coastal Boat Tensor Tactical 755.48 m 

4 Coastal Boat Tensor Navigation 35.94 m 

5 Coastal Aerial Vector Tactical 254.96 m 

6 Coastal Aerial Vector Navigation 103.43 m 

7 Coastal Aerial Tensor Tactical Filter Diverged 
8 Coastal Aerial Tensor Navigation 61.07 m 

9 Continental Aerial | Vector Tactical 514.53 m 

10 Continental Aerial | Vector Navigation 141.17 m 

11 Continental Aerial | Tensor Tactical Filter Diverged 
12 Continental Aerial | Tensor Navigation 205.95 m 

13 Deep-Sea Boat Vector 10X 194.90 m 

14 Deep-Sea Boat Tensor 10X 185.52 m 

15 Global Model Boat | Vector 10X 219.28 m 

16 Global Model Boat | Tensor 10X 38.10 m 
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V. Conclusion 


Through simulation, this research was able to prove the viability of using vector 
or tensor measurements of the magnetic anomaly field in a navigation system as 
an alternative to GPS navigation. Several variables were evaluated in the sixteen 


simulation cases and findings for these are listed below. 


Trajectory Position and Velocity. 


The velocity of the trajectory did not have a great affect on filter performance 
during simulations. This may be due to the fact that more signals are being brought 
into filter with vector and tensor measurements, so despite the measurements coming 
in more slowly, the higher number of signals allows the filter to accurately resolve a 
position and attitude solution. The navigation system was able to obtain navigation 
accuracies of 35.94 m over a one hour low-velocity boat trajectory and 61.07 m over 


a one hour high-velocity aerial trajectory. 


Measurement Types. 


The simulations demonstrate that navigation accuracy increases with an increas- 
ing number of measurements coming into the filter. As long as the filter was paired 
with a nav-grade INS, the navigation filter when using tensor measurements was able 
to achieve the lowest DRMS error (specifically in the coastal boat TEN-NAV case) 
and was able to consistently perform better than when using vector measurements. 
These are promising results, because the use of tensor measurements essentially can- 
cels out the effects of temporal variations in the measurements, which in turn reduces 
measurement errors in the navigation system. The ability to minimize measurement 
errors in å navigation system is critical to minimizing error in the position and atti- 


tude solutions. 
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The best case tensor filter results only occurred in combination with a navigation- 
grade INS. The tensor measurement filter was unstable in combination with a tactical- 
grade INS. Removing the temporal variation states from the model may actually have 
negatively impacted the filter. The filter was not expecting any bias errors because 
of the removal of these states; However, there may have been bias errors in the 
measurements due to other effects. This led to overconfidence in the filter estimates 


and in turn, filter divergence. 


INS Quality. 


The quality of INS used in simulations had a large impact on the overall filter 
performance. Generally, pairing the filters with a tactical-grade INS resulted in di- 
vergent runs. The filter performed best when paired with a navigation or 10X-grade 
INS. When paired with the 10X-grade INS, the navigation system was able to ob- 
tain navigation accuracies of 185 m DRMS over a 25 hour trajectory. This specific 
trajectory was over the deep ocean and used a map with lower frequency content. 
The TEN-10X filter was also able to obtain navigation accuracies of 38.10 m when 
paired with EMM global maps, which are of even lower frequency than the deep-sea 
NAMAD maps. We would expect the lack of signal in the map as well as lack of signal 
due to the low velocity of the ship to cause problems for the filter. The combination 
of the 10X-grade INS and the tensor measurements within the navigation system was 


able to overcome this challenge. 


Map Resolution. 


The simulations demonstrated that both navigation systems were able to perform 
with lower-resolution global maps from the EMM. When using either navigation sys- 


tem presented in this research, high navigation accuracy was not limited to navigation 
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using high resolution maps such as the NAMAD. For military applications, the results 
of this research open up the possibility of using magnetic navigation as an alternative 
to GPS navigation in any area of operations around the globe to estimate vehicle 


position and orientation accurately. 


Future Work. 


One aspect to consider for future work following this thesis, is running simulations 
on the navigation system with an alternate dynamics model. The EKF could possibly 
be improved when paired with the tactical grade INS and tensor measurements by 
introducing the error bias states back into the dynamics model. While the temporal 
variations are nearly canceled out in the tensor measurements, there could be more 
error from the effects of the vehicle than is modeled. More accurately modeling this 
error could improve filter performance, as the filter would place å more accurate 
amount of confidence in the measurements as opposed to the over-confidence that 
was displayed. 

Another option for dealing with this filter instability would be to re-design the 
filter to use an MPF as opposed to the EKF to better deal with the non-linearity of 
the measurement function. This would require a large amount of processing power 
given the computational intensity of dealing with the particles, but also computing 
the measurement Jacobian for each particle at each time step. This filter would be 
expected to have the highest accuracy if effort is put forth into minimizing processing 
requirements. 

Another consideration for future work following this thesis would be to move from 
simulation to real trials to assess the navigation system performance. While the 
trajectories used in the simulations were realistic, the performance of the navigation 


system using real measurements, maps, and trajectories could vary greatly depending 
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on the true accuracy of the measurements and maps. Vector magnetometers are a 
mature technology, however, they have not been widely used to conduct surveys 
and gather accurate vector magnetic anomaly maps. As these vector maps are not 
readily available, effort would go into obtaining accurate vector maps for the area 
traversed during real vehicle tests. The same is true for tensor maps however, the 
tensor measurement configuration is not as mature as the vector magnetometer and 
would require extensive calibration to ensure accurate tensor measurements. Once 
calibration is acceptable, the focus would be on collecting data to generate small 
tensor map tiles for real vehicle tests. Running trials with accurate maps and true 
measurements is the next step to determining if vector or tensor magnetic anomaly 


field measurements are å promising option as å GPS signal alternative. 
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